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Abstract. In the dynamic analysis of structural engineering systems, it is 
common practice to introduce damping models to reproduce experimentally 
observed features. These models, for instance Rayleigh damping, account for 
the damping sources in the system altogether and often lack physical basis. We 
report on an alternative path for reproducing damping coming from material 
nonlinear response through the consideration of the heterogeneous character 
of material mechanical properties. The parameterization of that heterogeneity 
is performed through a stochastic model. It is shown that such a variability 
creates the patterns in the concrete cyclic response that are classically regarded 
as source of damping. 

Keywords: damping ; concrete ; nonlinear constitutive relation ; material hetero- 
geneity ; stochastic field 


1. Introduction 

In the last few decades, a great deal of attention was paid to the comprehension 
and modeling of damping mechanisms in inelastic time-history analyses (ITHA) of 
concrete and reinforced concrete (RC) structures O section 2.4]. Figure[l| adapted 
from |33| . shows the uniaxial cyclic compressive strain-stress (E-E) response mea¬ 
sured on a concrete test specimen. Throughout this paper, the term “uniaxial” 
implies that there is only one loading direction and that the stress, respectively 
strain, of interest is the normal component of the stress, respectively strain, vector 
in the loading direction. In other words, when it comes to constitutive relation be¬ 
tween stress and strain, the work presented thereafter is developed in a ID setting. 
In figure the so-called backbone curve, which is the envelope of the response 
(dashed line), shows the following phases: (i) an inelastic phase with positive slope 
{E < 2.7 X 10“^ for that particular example, where E is the measured strain), 
and (ii) an inelastic phase with negative slope before the specimen collapses. For 
concrete, no elastic phase can really be identified, and hysteresis loops appear in 
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unloading-reloading cycles even for limited strain amplitudes. Other salient fea¬ 
tures include: (i) a residual deformation after unloading, and (ii) a progressive 
degradation of the stiffness (slope of the unloading-loading segments). The hys¬ 
teresis loops are one of the sources of the damping that is observed in free vibration 
recordings of concrete beams. Other sources include friction at joints m or at the 
concrete-steel interface in reinforced concrete [M] . These other sources of damping 
will not be discussed in this paper, where we will concentrate on material damping. 
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Figure 1. Strain-stress concrete experimental response in pseudo¬ 
static cyclic uniaxial compressive loading (adapted from [ 33 ] )• 5] 
and E are the homogeneous compression stress and strain in the 
concrete test specimen that are measured in the loading direction. 

S and E are spatial mean quantities in the sense that S is com¬ 
puted as the load in the hydraulic cylinder of the testing machine 
divided by the area of the specimen cross section, and E is com¬ 
puted as the displacement of the cylinder divided by the length of 
the concrete specimen. 

Most classical uniaxial constitutive models of concrete for numerical simulation 
do not dissipate any energy in unloading-reloading cycles (see e.g. figure[^[top left]). 
It is then common practice to add a viscous damping model (Rayleigh damping) 
to the inelastic structural model, to reproduce phenomena that are experimentally 
observed at the structural level (decreasing amplitude of displacements in free vi¬ 
bration). However, Rayleigh damping is well known to lack physical justification, 
even when care is taken to avoid generating spurious damping forces IZKIH]. 

Another class of approach aims at reproducing more precisely the features of 
Figure through elaborate inelastic constitutive relations. Figure shows typi¬ 
cal examples of uniaxial relations found in the literature. The relation described 
in |3| [top left] defines different response phases for different strain intensities, with 
an additional coefficient to control the loss of stiffness. The constitutive relation 
described in [23] [bottom right] comes from a formulation developed in the frame¬ 
work of thermodynamics with internal variables. It reproduces damping features 
reasonably well, in particular for higher amplitudes, but requires the identifica¬ 
tion of a rather large number of parameters. These first two types of relations 
are somehow defined by parts for different loading regimes. They hence require a 
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wider set of parameters and seem to contradict the seemingly smooth transition be¬ 
tween regimes observed experimentally. The constitutive relation described in |48j 
[top right] is heuristically defined from a database of experiments. It reproduces 
unloading-reloading hysteresis mechanisms, but lacks a theoretical basis. Finally, 
the relation described in [32] [bottom left] is based on a physical model of damage 
and friction. It manages to dissipate energy in unloading-reloading cycles, but the 
lack of obvious physical meaning for some parameters can render their identification 
difficult. 



Figure 2. Typical strain-stress relations in pseudo-static cyclic 
compressive loading for different models: |9] [top left], [48] [top 
right], [32] [bottom left], and [23] [bottom right]. 

The main purpose of this paper is to present a multi-scale stochastic nonlinear 
concrete model that can be accommodated in an efficient structural frame element 
(fiber element), and that participates to the overall structural damping in dynamic 
loading. In particular, this implies the developed concrete model be capable of rep¬ 
resenting hysteresis loops in unloading-reloading cycles at macro-scale (the scale 
where such behavior as in Figurecan be observed). In this work, this is achieved 
by the introduction, at an underlying meso-scale, of spatial variability in the pa¬ 
rameters. At meso-scale, an elasto-plastic response with linear kinematic hardening 
and heterogeneous yield stress is considered. This choice is mainly driven by its 
simplicity and its relevance is illustrated in the numerical applications. 

The main issue with modeling the heterogeneity of the yield stress lies in the pa¬ 
rameterization. On the one hand, local information on the heterogeneity of concrete 
is available at a scale that we wish to avoid (because of the associated computa¬ 
tional costs). On the other hand, identification becomes extremely difficult when 
very fine models are considered. We therefore choose to model the heterogeneity of 
the yield stress through a stochastic model. Hence, only three parameters control 
that heterogeneity: a mean value, a variance, and a correlation length. The choice of 
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parameterizing the fluctuating field of constitutive parameters by statistical quan¬ 
tities means that there might be fluctuations in the quantities of interest measured 
for different realizations of the random model. However, as will become apparent 
in the examples, some sort of homogenization comes in and these fluctuations can 
rightfully be ignored. 

Several authors in the literature have considered random models of fluctuating 
nonlinear materials [111 HOJ [H US], in particular for concrete [HI [331 HSl ESI ESI 
[SOI [30] or in the context of dynamic analysis [311 UHl ES]- We consider here a 
modeling framework that is a combination of ingredients found in several previous 
papers |23|S1|3I1|30|, with a fluctuating yield stress modeled as a random field with 
non-zero correlation length. However, the objective in these papers was to assess 
the influence of parameter uncertainty on some quantity of interest. An objective 
with the current paper is to observe the effect of randomness at a meso-scale on the 
nonlinear stress-strain relation at macro-scale. The work herein presented should 
therefore be seen as an innovative proposal for parameterization of a nonlinear 
stress-strain relation. 

In Section we recall the theoretical formulation of the inelastic beam model 
that will be used throughout this paper. The stochastic multi-scale constitutive re¬ 
lation developed to represent concrete cyclic behavior in reinforced concrete frame 
elements is introduced in section [Sj Concrete behavior at macro-scale is retrieved 
from the description of a meso-scale where elasto-plastic response with linear kine¬ 
matic hardening and spatially variable yield stress is assumed. In particular, we 
emphasize in sections |3.3| and [3.4| the heterogeneity of the yield stress and the pa¬ 
rameterization of that heterogeneity through a random model. In section we 
report on the limiting case of vanishing correlation length and monotonic loading, 
for which several results can be derived analytically. Section presents numerical 
applications of the model in the context of dynamic structural analysis of reinforced 
concrete frame elements. 


2. 2D CONTINUUM Euler-Bernoulli inelastic beam 


Classical displacement-based formulation has been retained here although other 
mixed formulations can in certain cases show better performances |46j . For the 
sake of conciseness, we present the beam element in the 2D case, extension to 3D 
is straightforward. 


2.1. Euler-Bernoulli kinematics. We define the continuum beam H = {x G 
G [0, L]; X2 G [—h/2, h/2]; 0:3 G [—w/2,w/2]}. Such a beam has length L 
and uniform rectangular cross-section S of size w x h. We consider an orthonormal 
basis (ii, i 2 , is) of so that any material point in space x = X]i=i the 2D 

setting adopted here, Euler-Bernoulli kinematics can be written at any point x G H 
and at any time t G [0, T] as 


( 1 ) 


u(x,t) 


Ui (x, t) = uf (xi ,t) - X20§{xi, t) 
U2(x, t) = U2{xi,t) 


with 


0iixi,t) = 


du2{xi,t) 

dxi 


Ui and U 2 are the longitudinal and transversal components of the displacement 
vector u at any point in the beam, uf and uf are rigid body translations and Of 
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is rigid body rotation of section S at position xi along the beam axis. Thus, uf, 
rtf and 9§ only depend on a;i. 

For small transformations, strain tensor reads E = ^ (D(u) + D^(u)), where 
D(-) = Ell & ® b, with ® the tensor product and the transpose operation. 
Then, defining the axial strain = duf/dxi and the curvature = d'^u^ldxj, 
it comes: 

(2) E(x, f) = i:;(x,t) ii (g) ii 


with 

(3) E{-x,t) = e^(xi,t) - X 2 X^ixi,t) . 


2.2. Variational formulation. Suppose B is loaded with forces per unit length 
of beam b and concentrated forces applied at beam ends. A variational form of the 
problem of hnding u such that beam equilibrium is satisfied is: find n such that 

(4) 0 = y |y (^e‘^ — X2Sx^) E dxi — J • b dxi — , 

where 5n is any kinematically admissible displacement field E = S • ii (g) ii with S 
the stress tensor, • a matrix product here, and 6Ilbc the potential for the forces at 
beam ends. 

Introducing the normal and bending forces in the beam cross-sections as 


( 5 ) 


N= E diS and M = — X 2 'S, dS , 


equation Q can then be rewritten as 


( 6 ) 


0 = / Se^ ■ q dxi — / i5u • b dxi — , 


where q = {N,M)'^ and e*^ = 


2.3. Inelastic constitutive behavior. Cross-section inelastic constitutive response 
q(e‘^; t) is thereafter represented using uniaxial material constitutive response E(E; x, t) 
integrated over the cross-section, rather than using a direct relation between section 
displacements and forces. This approach leads to what is often referred to as fiber 
beam element. 

With A denoting an increment of some quantity, we introduce the tangent mod¬ 
ulus D as AE = D X AE, that is, from equation ([^, 

(7) AE,{x,t) = D{x,t) {Ae^{xi,t) - X 2 Ax^{xi,t)) . 

Introducing relation Q in ([^, beam section inelastic constitutive equation reads 
Aq = K'^Ae'^ with tangent stiffness matrix 


L./5 


I 

-X2 


D{x,t) (1 


—X 2 ) dS 


( 8 ) 


K'5(xi,t) 
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2.4. Numerical implementation (structural level). The finite element method 
is used to approximate the displacement fields. Classically, we have for each ele¬ 
ment u(x, t) = N(x)d(t) and e‘^(x, t) = B(x)d(t), where the vector d gathers the 
displacements wf,uf ,03 at the element nodes, and matrices N and B gather the 
classical shape functions for Euler-Bernoulli kinematics. Choosing, in equation Q, 
(5u = Ni5d and Se'^ = B(5d, and then linearizing the resulting relation, we have 

(9) K«Ad«=r«, 

where K = /^ B^K'^Bdxi and r = f — B^ qdxi- f is the vector of nodal 
forces calculated from b and from any concentrated force applied at a beam ends. 
Subscript n refers to any time step in the loading history; superscript k refers to 
any Newton-Raphson iteration. 

For any integrable function g, fine integrals are numerically approximated as 
J^g{x)dx « gixi)Wi where subscript I refers to a quadrature point and 

Wi denotes quadrature weight and length. Section integrals are estimated as 
Is, ffMdSi « X;f=i ), where is the section area of the so-called fiber 

F and xf is the position of the fiber centroid in the control section Si at quadrature 
point 1. 


3. Multi-scale uniaxial cyclic model for concrete 

In this section, we present the core objective of the paper, which is a nonlin¬ 
ear uniaxial constitutive model capable of representing salient features of concrete 
compressive response in cyclic loading. It is based on a simple local (meso-scale) 
elasto-plastic constitutive relation, for which the yield stress is modeled as a random 
field. The spatial fluctuations of the yield stress induce at macro-scale constitutive 
relation E(i?) which resembles that encountered experimentally. Again, as already 
stated in the introduction, the idea of considering an elasto-plastic constitutive re¬ 
lation with fluctuating yield stress is not novel per se [IS [3 mi |3Q]. It has been 
proposed to assess effects of uncertain parameters on model outputs of interest for 
engineering practice while our aim here is to stress on the fact that this can be seen 
as a way of parameterizing material nonlinear constitutive relations. In particu¬ 
lar, it will be shown that, in some circumstances, even if randomness is present at 
meso-scale, our model can predict non-random outputs. 

3.1. Meso- and macro-scale modeling of concrete. Concrete is a heteroge¬ 
neous material (see figure]^. Two scales are classically considered for its modeling: 
(i) a micro-scale at which each phase (aggregate, concrete, cement paste) is clearly 
identified and modeled with its own constitutive relation; and (ii) a macro-scale at 
which concrete is considered as homogeneous. The macro-scale is the relevant scale 
for structural engineering applications but the behavior at that scale is strongly 
influenced by phenomena occurring at the micro-scale. In particular, the geome¬ 
try of the phases is important as it controls in a large part local concentrations 
of stresses. As pointed out in the introduction to this paper, formulating and im¬ 
plementing concrete constitutive laws at the macro-scale can then turn out to be 
challenging, even in the uniaxial case, and especially when it comes to accounting 
for material energy dissipation sources. We follow here another path, considering 
a meso-scale at which the parameters of the constitutive relation are assumed to 
vary continuously. This scale is intermediary between the macro-scale, at which 
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the parameters are homogeneous, and the micro-scale, at which the parameters are 
discontinuous. 



Figure 3. Polished concrete section where the two phases human 
eyes can see are represented: aggregates (crushed gravel and sand) 
and cement paste in-between. 


For practical implementation, the heterogeneity will be conveyed in our model 
by the fluctuations of a random field p(x, d), where 9 represents randomness. Con¬ 
sistently with the fiber beam formulation presented in the previous section, we 
consider a mesh of fibers F spanning beam cross-sections S. These fibers have a 
centroid located at position and a cross-section denoted by TZ. In the spirit of 
strain-controlled tests on concrete specimens (see figurealong with its caption), 
strain E is assumed homogeneous over TZ: 

(10) e{x,t) = E{x[,t) Vx G 7?, , 


and stress at macro-scale S is computed as the spatial mean stress over TZ: 

( 11 ) E{E;9) = ^J^a{e-,e)dTZ, 


where \TZ\ denotes the area of the fiber section TZ, a and e the stress and strain at 
meso-scale in the uniaxial case. That is, analogously to E and E, e = e • ii ® ii and 
(T = cr • ii 0 ii, where e and rr are the strain and stress tensors at meso-scale such 
that the constitutive model is considered in a ID setting. Besides, we enhance the 
fact that E is computed as the spatial mean of cr(x) and not as the sample mean of 


(t( 0). Local constitutive relation presented in section 3.2 below governs the relation 
between heterogeneous stress field a{x,t) and strain field e(x, t). 

Then, we introduce the tangent modulus at meso-scale ® as Act = D x Ae. 
According to equations 0, ([^ and 0. we have the tangent modulus at macro¬ 
scale 


( 12 ) 


D(x, t; 6») = -j^ ^ £>(x, t; 9)dTZ 


We will see in the examples below that, for a wide range of relative correlation 
length and size of the section TZ, E and D do not depend on 9. In that case. 
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even though the meso-scale model of concrete is stochastic, the resulting macro¬ 
scale model is deterministic, and independent of the actual realization of the local 
parameters that is being considered. 

Finally, we recall that concrete specimens exhibit quasi-brittle behavior in ten¬ 
sion with tensile strength generally 10 times smaller than compressive strength. 
Besides, we point out here that it is at macro-scale that these notions of compres¬ 
sion (S < 0) and tension (S > 0) are relevant for the modeling presented in this 
work. 


3.2. Model of the uniaxial cyclic behavior at meso-scale. We now concen¬ 
trate on the local uniaxial cyclic inelastic constitutive relation that will be consid¬ 
ered in this paper at meso-scale. We define it to be a simple elasto-plastic model 
with linear kinematic hardening, as illustrated in Figure We provide here, in 
the setting of computational inelasticity mmn, the assumptions and resulting 
equations corresponding to this relation: 

(i) The total deformation e is split into elastic (e®) and plastic (e^) parts: 

(13) e = + eP . 

(ii) The following state equation holds (upper dot denotes derivative with re¬ 
spect to time): 

(14) d = Ce® , 


where C is the elastic modulus. 

(iii) We impose that the stress a corrected by a, the so-called back stress due 
to kinematic hardening, satisfies yielding criterion 


(15) 


(j)P = \a + a\ — ay < 0 , 


where Uy > 0 is the yield stress. As yielding function cj)P{a,a) is negative, 
the material is elastic; otherwise, plasticity is activated and the material 
state evolves such that the condition (j>'^{a, a) = 0 is satisfied. 

(iv) A change in can only take place ii (j)^ = 0 and yielding occurs in the 
direction of cr -I- a, with a constant rate 7 ^ > 0: 


(16) 


J 7 Psign(CT -I- a) 

I 0 


if (j)^{a, a) = 0 
otherwise 


jP is the so-called plastic multiplier. 

(v) With H the kinematic hardening modulus, the evolution of a is defined as: 


(17) 


a = —He^ = —7^7Jsign((T -I- a) . 


Accordingly, admissible stresses cr and a remain in the set /C® = {(cr, a)\ (jiP < 0} 
and two kinds of evolutions are possible: 

(i) If (ct, a) G JC^ = {{a, a) | (()^ < 0}, the response is elastic: 


( 18 ) 


eP = 0 


= Ce . 


(ii) If (cr, a) G i9/C® = {(cr, a) \ cj)P = 0}, any evolution is possible only if = 0: 

CH . 


dr dr. 


C'sign(cr-I-a)e 
^ C+H C+H 
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Figure 4. Compressive cyclic behavior at meso-scale. The yield 
stress cry(x) fluctuates over the concrete section TZ: two local re¬ 
sponses at two distinct material points Xi [left] and X 2 [right] are 
represented in the figure. During elastic loading/unloading, the 
slope is C; in yielding phases, the slope is < C. With this 
type of model first yielding occurs once a = cfy and the elastic do¬ 
main keeps constant amplitude 2ay. For low values of (Ty, yielding 
can be observed in compression both during loading and unloading 
[right]. 


It is then possible, from equations (18) and (19), to give the expression of the 
tangent modulus S: 


( 20 ) 


a = 2)e with 2) = 


C 

CH 

C+H 


if {a, a) G JC^ 
if (a, a) G 


It should be reminded at this point that, as mentioned in the introduction and 
illustrated in Figure the yield stress is assumed heterogeneous. The relation 
presented here is therefore defined between stress and strain in each point in space 
with a different yield stress. 


3.3. Description of the yield stress random field. In this section, we describe 
the choice that is made for the modeling of the heterogeneous yield stress: the yield 
stress is represented by a 2D log-normal homogeneous random field over the con¬ 
crete area TZ. We note here that, to the best of our knowledge, there currently exists 
no experimental dataset of local stress-strain uniaxial concrete responses recorded 
at many points over a concrete area. Here, the choice of using random fields to con¬ 
vey heterogeneity of the yield stress is mainly motivated by the effectiveness of the 
method. We hope that this proposed interpretation of concrete meso-structure will 
foster interaction between numerical and material scientists and help designing ex¬ 
perimental investigations that would eventually support or invalidate the numerical 
model we propose in this paper. 

Let us then consider a probability space (0,D,Pr), where D is a tr-algebra of 
elements of 0 and Pr is a probability measure. The 2D random field of yield stress 
is constructed as a nonlinear point-wise transformation [13 = f(G(x;6l)) 

of a homogeneous unit centered Gaussian random field G(x; 6) with given power 
spectral density (PSD) Sccif^)- The PSD is chosen here as the product of triangle 
functions with identical properties in the two orthogonal directions of the 2D plane. 
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denoted by the subscript 1 and 2 throughout sections |3.3| and |3.4| 


( 21 ) 


'S'gg(«^) 



where A(k) = 1 — |k| if |«;| < 1 and cut-off wave numbers = Ku ,2 = I’^u- For 
wave numbers above the cut-off the spectral density vanishes. In the spatial 
domain, this PSD corresponds to the following autocorrelation function (the Fourier 
transform of S'gg(<«)): 

(22) i?GG(C) = sW (^Ci) sW (^C2) , 

where sinc(a;) = sin(7ra;)/(7rai). The random field G(x;0) therefore fluctuates over 
typical lengths £c,i = ^c ,2 = the so-called correlation length. 

The nonlinear point-wise transformation f controls the first-order marginal dis¬ 
tribution of ay. In particular, it controls the desired expectation m and variance 
of the yield stress homogeneous random field. In this paper, we choose to consider 
a log-normal first-order marginal density, to ensure that the realizations are almost- 
surely and almost everywhere positive, as expected. The nonlinear transformation 
is then given by: 


(23) 
where 

(24) 


6y(x, 9) = exp(mG + sg x G(x, 6*)) > 0 , 



and 



Other first-order marginal densities could be considered, for example using the 
maximum entropy principle [m Ea Ea sni HI] or Bayesian identification 
Also, the PSD function is translated by the nonlinear transformation f so that the 
PSD of the yield stress and of the underlying Gaussian field are different, with 
possible incompatibilities with the chosen first-order marginal density HZlISIlEHl- 
These important but technical issues go beyond the scope of this paper and will 
not be further discussed here. 



Figure 5 . Realizations of log-normal random fields over a square 
of size d for different correlation lengths: [left] = 1, [center] 
idd = 0.1, [right] G/d = 0 (white noise). Coordinates xi and 
xs, previously introduced in the description of the beam element 
in section are reused here to recall that the random fields are 
generated to parameter heterogeneous yield stress over beam cross- 
section areas. 
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3.4. Numerical implementation (at each quadrature point in each beam 
fiber). At each material point xf (quadrature point I, beam fiber F), relations (101, 
0 and ( [T^ have to be calculated. 

On the one hand, Gaussian random field G is digitized using the spectral rep¬ 
resentation method, in its FFT implementation [391 na 1^ . As an illustration, 
considering a 2D random field with identical properties in orthogonal directions 1 
and 2 (see [ID] for more details): 


M-l M-1 

(25) G(piAa:,p 2 Aa:; 6 ») = Re ^ ^ ( 6 ») exp (^ 2 f 7 r ) 

ni=0 712—0 

+R„,„,( 0 ) exp (2*71 (^-^))) 

where = - 1 , (pi,_P 2 ) G [ 0 ,... ,M - 1 ]^, 

Bnin2 = 2AKy SoainiAK, n2 Ak) exp(*(/)„j„2 (6»)) 

(26) ScainiAK, -n2AK) expii'ipriinA^)) 


and Ak = k„/A (A S N* —>■ oo), M > 2A, are independent 

random phase angles uniformly distributed in [0,27r]. The resulting 2D random 
field is periodic with a two-dimensional period Lq x Lq with Lq = MAx = 2 tt j 
Ak. Realizations of log-normal random fields with different correlation lengths are 
shown in figure 

On the other hand, a generic concrete section TZ is built as a square with edge 
of length d and TZ is meshed by a square grid of Nj identical squares. Then, the 

mesh size is d/Nj and, for any integrable function g, g{x)dTZ « ^ X(/=i 

where x-f is the position of the centroid of the f-th mesh over TZ. 

Then, digitized random field Gy is mapped onto the x-^’s over TZ. To this purpose, 
we impose Lq > d, that is \TZ\ is smaller or equal to a period of the random field, and 
mapping is performed according to the following method. First, Nf is calculated 
as: 


(27) 


d = Int 



Ax -I- Res => Nf = 


lnt{d/Ax) if Res = 0 

Int(d/Ax) -I- 1 otherwise 


Then, at the A| points x-t S 77, &y{xd]9) is calculated as the linear interpolation 
of the four digitized values of 6y(x;0) in Jx-t — Ax,x^ + Ax]^, as illustrated in 
figure 

With the spatially variable yield stress now known at each point x-f, / G [1,.., A|] 
in TZ, the equations presented in section can be solved to update the variables at 
meso-scale. This is done numerically at each of the A| positions following classical 
return-mapping computational procedure mnn]. 

Finally, as a transition from compression to tension is detected during global 
Newton-Raphson iterative process to solve structural equilibrium equations, that is 
> 0 while < 0, a local Newton-Raphson process is implemented to find 
the strain for which = 0, to update the meso-structure accordingly, 

and to set = 0 and = 0. 

Before observing on numerical tests the shape of the stress-strain curves obtained 
with this model, we turn to the simple case of vanishing correlation length (f^ —t 
0). The interest of this particular case is that some analytical expressions can be 
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Figure 6. On the one hand, the random field is digitized on a 
square grid (dashed lines) of (M + 1)^ points — here M = 6 — 
spanning an area of size LqX Lq, with Lq = M x Ax. On the other 
hand, generic concrete square section TZ has an area of size d x d 
{d < Lq) that is divided into Nj identical meshes (plain lines) — 
here Nf = 5. The values of the random field at the positions x-^ 
occupied by the centroids of the Nj meshes (x) is calculated as 
the linear interpolation of the four surrounding digitized values of 
the random field (o). 


derived, so that discussion is more straightforward. The more general case with 
finite correlation length will be considered later in Section [5T| 


4. A PARTICULAR CASE: VANISHING CORRELATION LENGTH AND MONOTONIC 

LOADING 

4.1. Preliminaries. The case of vanishing correlation length along with uniaxial 
cyclic loading has been treated in a general setting. Indeed in [21], the stress- 
strain uniaxial response is given as a probability density function (pdf) of stress 
with respect to the time-dependent strain and a second-order exact expression of 
the pdf evolution is computed solving the Fokker-Planck-Kolmogorov equation that 
governs the problem. This latter method is valid for monotonic as well as cyclic 
loading. Hereafter, the validity of the results is limited to monotonic loading, but 
the problem is cast in a different and simpler mathematical setting that can be 
solved analytically. These analytical developments shed light on some capabilities 
of the model introdnced in the previous section and that will be retrieved in the 
more general case of non-zero correlation in Section |5.1[ 


4.2. Constitutive response at macro-scale. Let respectively denote 7^® and 
TZP the shares of a fiber cross-section that remain elastic and yield. According 

TZ^{t;9) = {x G 7^ I X>(x, t;0) = C} and 
D(x,t;6») = ^/{C+H)}. Also, 7^®^7^P = 0 and 7^ = 7^®U7^P. 


to the developments in section 3.2 
7^P(^;6») = {x G 7^ 


Note that 7?.® and TIP are time-dependent because S) depends on the loading history. 
We denote by | • | the area of •. Then, considering a subset A of TZ, we have, Vx G TZ, 
the probability measure Pr[x G A] = |A|/|7^|. 






























13 


Using the fact that \R.\ = 
macro-scale in equation (121 as: 


(28) 




\'RP\, we first rewrite the tangent modulus at 
C 


H 


c + H\ \n\ 


C + H 


We now seek an explicit expression for |7?.®|/|7?.|. 

First, suppose the state of the material is known at time tg, then we define the 
trial stresses 


(29) 


cr‘’’(x,t) = (To(x) + C(e(t) — eg) and a*^{x,t) = ao(x) ; 


where subscript 0 refers to time tg. In the particular case of monotonic loading, a 
necessary and sufficient condition for x to be in 7?,^ at time t > to is (^^’*’'(x, t) > 0, 


that is a-y(x) < |(t*’'(x, t) -|- ao(x)| (see equation 15). We then have: 


(30) 


|7?.®|/|7?.| = Pr[x e TV'] = 1 - Pr[6j^(x) < |cr‘’'(x, t) + ao(> 


Then, in the particular case of vanishing correlation length, the random variables 
&y{x) are independent and identically distributed over TZ. For the log-normal 
distribution assumption made throughout this work, it means that the cumulative 
density function of &y{x) is, Vx G TZ: 

(31) -^eUxlK) =Pr[6^(x) < ay] = ^ (^1-f erf > 

where erf is the so-called error function. 

Finally, for the sake of simplicity and without any loss of generality, we assume 


CTo = ao = Co = 0. Accordingly, and using equations (29), along with (24) to 


replace mo and sq by the mean m and standard deviation s of the homogeneous 
log-normal random field &y, it comes: 


(32) 


^ = l-^eUx)(|C^eWI) = ^ 


V 


1 — erf 




V 






Equations { 28) and (32) are used to plot figurej^where the response of the model 
at macro-scale is shown for different sets of mean and variance parameters for the 
log-normal random yield stress field &y. 

4.3. Asymptotic response of the model at macro-scale. The following as¬ 
ymptotic behaviors can be observed at macro-scale: 


(i) Suppose jrr? approaches 0. Then, according to equation (32), |7?.®|/|7^| 
approaches the Fleaviside’s function 'H{m— ]CE{t)]), that is |7?.®|/|7^| = 0 if 
]CE{t)] > m and |7^®|/|7?.| = 1 if ]CE{t)] < m. According to equation (28), 
the model response at macro-scale is then as follows: 


(33) S(t) = D(t)E{t) where D = 

(ii) If ^jrr? —>■ oo, then |7^®|/|7^| 
H)E{t). 


if |C'F;(t)|<TO 
if ]CElt)]>m 


C 

CH 
C+H 

> 0 and consequently E(t) 


CH/{C ■ 
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Figure 7. Monotonic macro-scale response of the model with 
vanishing correlation length. Parameters C = 30 GPa and 
H = 10 GPa are used. [top] m = 30 MPa and s/ 
m = 10“"^, 0.4,2,10,10®; [bottom] s = 30 MPa and s/m = 
10-3,0.5,1,2.5,103. Plain curves show asymptotic behaviors as 
s/m —0 or - 1 - 00 . 

(hi) Now with finite and non-zero s^/m^: 

(34) ±{t) = D{t)Eit) where ' 

These asymptotic responses at macro-scale are illustrated in figure(plain lines). 

5. Numerical applications 

5.1. Concrete uniaxial compressive cyclic response at macro-scale. First 
numerical applications aim at demonstrating the capability of the model introduced 
above in section to represent the response of concrete in uniaxial compressive 
cyclic loading. Five model parameters need to be considered: elastic and harden¬ 
ing moduli C and H, along with mean m, standard deviation s and correlation 
length ic used to build realizations of a homogeneous log-normal random field that 
parameterizes the fluctuations of the yield stress ay over beam sections. 

The effects of to, s and ic on the material response at macro-scale will be further 
investigated below. Right now however, we set: 

• C = 27.5 GPa, which corresponds to the elastic modulus measured on 
specimens made of the concrete actually cast to build the frame element 
used in the next numerical application (Section |5.2[ ). 

• H = 0 according to both (i) the fact that H controls tangent modulus at 
macro scale as strain becomes large (see equation @)) , and (ii) that we 
seek a numerical response that ultimately exhibits null tangent modulus 
in monotonic loading at macro-scale. We anticipate here stressing that 
the model developed in previous sections is not capable of representing 
the softening phase as strain increases while stress decreases (non-positive 
tangent modulus). 
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Figure 8. Sample mean (thick plain line) plus/minus standard 
deviation (boundaries of the shaded areas) monotonic response at 
macro-scale computed from a sample of 100 realizations of the 
material structure at meso-scale with to = 30 MPa and s/to = 1 
for the log-normal marginal law. Meso-structures are generated 
with different correlation lengths: [left] Icjd = 0.1, [center] l^j 
d = 0.2, [right] £c/d = 0.4. Cyclic response for one particular 
realization of the meso-structure is also shown (thin plain line). 


5.1.1. Influence of ic on the macroscopic response. We first illustrate how £c in¬ 
fluences the macroscopic response by considering the three following cases: (i) id 
d = 0.1, (ii) £c/d = 0.2 and (iii) £c/d = 0.4. For each of these three cases, we 
take = 10 and M = 64, that is Ax/d = 0.016 and Nf = 64 (see equation (27l). 
We recall that the random field characteristics are taken as identical in both space 
directions {d = £c,i = £c, 2 , N = Ni = IV 2 , ■ ■ ■ )• Besides, a sample of 100 indepen¬ 
dent homogeneous log-normal random fields with targeted mean to = 30 MPa and 
coefhcient of variation s/to = 1 for the marginal log-normal law is generated for 
each case. 

Resulting material responses at macro-scale are shown in figure [^ A first ob¬ 
vious observation is that model response at macro-scale is much richer than at 
meso-scale (see figure |^. We can then observe that sample mean response (thick 
line) is not sensitive to the correlation length. However the variability of the macro¬ 
scopic response from one realization of the meso-structure to another depends on 
the correlation length: it is almost null for £c/d = 0.1 while it is enhanced as £c 
increases. The area TZ defined with £c/d = 0.1, = 10 and M = 64, is statistically 

representative in the sense that there is almost independence between the random 
realization of the meso-structure and the response at macro-scale. We remark here 
that the model is capable of representing variability from one concrete sample to 
another and that this variability can bring information on the correlations in the 
meso-structure. Suppose indeed that we had 100 concrete samples and a variability 
of the responses at macro-scale close to that shown by the grey area in figure [^ 
[center] for instance. Then, the meso-structure of the tested concrete would be best 
represented by the ratio £cld = 0.2. Consequently, the sample standard deviation 
can bring information on the actual correlation length. 

Finally, we can notice that the shape of the cyclic response (thin line) represents 
most of the salient features exhibited experimentally in uniaxial compression test 
for concrete (remember figure [^. We point out here that strength degradation 
(softening) along with stiffness degradation (observed experimentally one cycle after 
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another) are not represented by this model. However, a key point for representing 
material damping is the capability of the model to generate local hysteresis loops 
in unloading-loading cycles. 



Figure 9. Sample mean (plain line) plus/minus standard devia¬ 
tion (dashed lines) response at macro-scale computed from a sam¬ 
ple of 100 different realizations of the material structure at meso- 
scale. Meso-structures are generated with different targeted mean 
m and coefficients of variation s/m for the log-normal marginal 
law: [left] m = 30 MPa and s/m = 0.1,1,3; [center] s = 30 MPa 
and s/m = 0.6,1,3; [right] s/m = 1 and m = 10, 30, 50 MPa. 

5.1.2. Influence of m and s on the macroscopic response. We use here, beside N = 
10 and M = 64 (TV/ = 64), Ic/d = 0.1 so that little variability is expected to be 
observed in the model response at macro-scale from one realization of the meso- 
structure to another (see figure]^ [left]). Figurej^shows material response at macro¬ 
scale for different sets of mean m and standard deviation s of the log-normal random 
field that conveys spatial variability in the material structure at meso-scale. 

It can be observed that for a small value of s, response approaches bi-linear elasto- 
plastic behavior (actually perfectly plastic because H is set to zero here) where there 
is almost no hysteresis observed during unloading-loading cycle (figure [left] and 
[right]). This comes from the fact that, if s approaches 0, there is almost no spatial 
variability of the yield stress because it is almost homogeneous over TZ and takes 
values close to its mean m; then the response at macro-scale coincides with that 
at meso-scale (elasto-plasticity with H = 0 here). This is also in accordance with 
what was shown already in figure 

Responses shown in figure [left] lie in-between this latter extreme case and 
the other extreme case of s approaching infinity. In this situation, log-normal 
distribution approaches 0 all over the positive real semi-line and, consequently, 
plasticity is activated almost everywhere over TZ resulting in a macro-scale response 
that is perfectly plastic without elastic phase (that is here S = 0 for any E because 
H = 0). 

Also, it is shown in figure [^ that the value of the strain E at which stress S 
reaches zero when unloading (residual plastic deformation) is much more sensitive 
to parameter m than s. Furthermore, the thickness of the hysteresis loops obviously 
depends on the so-called coefficient of variation s/m but it is not clear whether it 
is more sensitive to either of the two parameters. Finally, the variability in the 
sample of responses at macro-scale increases with s/m and is more sensitive to m 
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than s, at least as far as the range of values chosen here for both parameters is 
concerned. 

5.2. Damping in a reinforced concrete column in free vibration. We now 

show how the material model developed in the previous sections can be used to 
represent the experimental backbone curve of a concrete test specimen in uniaxial 
loading. Then, we implement this material law in the fiber beam element presented 
in section and show how damping is generated in a reinforced concrete (RC) 
column in free vibration. The observed damping does not result from the addition 
of damping forces in the balance equation of the RC column but from the hysteresis 
loops in the concrete material law at macro-scale. 

5.2.1. Geometry of the column and loading. The column considered here corre¬ 
sponds to the Ist-floor external column of the ductile {R = 4) RC frame tested 
in nails]. The loading is however here different: it consists of a mass M = 500 kg 
imposed step by step and kept constant while the column oscillates in free vibration 
consequently to a horizontal force F{t). The geometrical and loading characteristics 
of the column are depicted in figure 




z 


Section area: 

^ 0.14x0.18m 

Steel bar area: 

/4s = 

Section A-A 



Figure 10. Geometry and loading of the column. 


5.2.2. Conerete eonstitutive model. In [25| . the monotonic uniaxial response (back¬ 
bone curve) of the concrete cast to build the RC column is detailed. It is used 
here as the baseline to identify the parameters of the concrete model. According 
to this report, we set C = 27.5 GPa, H = 0; then we use = 10, M = 64 
{Nf = 64) and ic/d = 0.1 to ensure a response at macro-scale that is almost in¬ 
dependent of the realization of the meso-structure; finally, m and s are identified. 
Figure 11 shows model response macro-scale (plain line) with m = 30.5 MPa and 
sjrfi = 0.943. Comparing the numerical backbone curve (solid line) and the experi¬ 
mental response (dashed line), this figure illustrates the capability of the developed 
numerical model to represent actual experimental concrete monotonic response, at 
least as far as the monotonic behavior is concerned. Note that no experimental 
data was available for the cyclic behavior. 


5.2.3. Steel cyclic model. Young modulus Cg = 224.6 GPa, yield stress Sy = 438 
MPa and ultimate stress = 601 MPa have been experimentally measured during 
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Figure 11. Sample mean (—) response at macro-scale obtained 
numerically from a sample of 2,000 meso-structures, along with 
backbone curve (- -) recorded during uniaxial test on a specimen of 
the concrete used to build the RC column of interest here. Targeted 
mean and standard deviation of the marginal log-normal law are 
m = 30.5 MPa and s/m = 0.943. No experimental data is available 
for the cyclic behavior. 


uniaxial tests on longitudinal steel rebars An elasto-plastic model with kine¬ 
matic hardening is used to represent steel response in cyclic loading. The model 
implemented with these latter measured parameters is shown in figure [T^ 



Figure 12. Numerical cyclic response of a steel longitudinal rebar 
used to build the frame. Cyclic behavior has not been observed 
experimentally. 


5.2.4. Free vibration - Structural damping. Those concrete and steel uniaxial con¬ 
stitutive models are implemented in the fiber frame element presented in section 
The column is modeled with one frame element with Ni = 2 control sections and 
Np = 6 fibers (actually layers here in the case of a 2D problem). As already men¬ 
tioned in section |5.2.1[ the mass M = 500 kg is imposed step by step and kept 
constant while the column oscillates in free vibration consequently to the horizon¬ 
tal force F{t). The column possibly exhibits nonlinear response while the mass M 
is applied and while F(t) increases from 0 to Fg. Figure 13 shows typical column 
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top-displacement time histories for two different values of Fq. It can be observed 
that damping depends on the amplitude of the oscillations: the column is clearly 
damped for Fq = 15 kN (grey curve) while damping is much lower for i^o = 5 kN 
(black curve). One can also notice the different vibration periods for both hori¬ 
zontal forces; this is due to the fact that the larger force activates some nonlinear 
mechanisms in the structure, which leads to an elongation of the structural vibra¬ 
tion period. We finally stress again here that there is no damping force added 
in the dynamic balance equations, such as for instance Rayleigh damping forces: 
the damping effect shown in figure only comes from the hysteresis loops in the 
concrete response during unloading-reloading cycles. 



time [s] 


Figure 13. Top displacement time history in free vibration for 
Fq = 5 kN (black) and Fq = 15 kN (grey). Mass M and horizon¬ 
tal forces Fq are applied step-by-step during the first and second 
seconds, then horizontal force F abruptly drops to zero and the 
column oscillates in free vibration. 


We now define what we will refer to as “viscous-like damping ratio” and hereafter 
denote by Considering the column top-displacement time history Xtop{t) in free 
vibration, we appeal to the so-called log-decrement method to evaluate the modal 
damping ratio (see e.g. n §4.6]): 


(35) 


1 ^ 
27riV, 


and are the amplitudes of any two peaks separated by = 

N 2 — Ni cycles. It is worth recalling here that this is only valid in case damping is 
linear viscous, which in our case is not necessarily the case. Indeed, equations (35) 
comes from the assumption that the envelope of the decaying top-displacement is 
described as Xtop{t) = XQe~^^'"f* with / the modal frequency. Hence the terms 
“viscous-like” to characterize the calculated damping ratios. 

Based on the top-displacement time histories in figure figure [Td] shows how 

decreases throughout free vibration time history for both values of Fq. Viscous-like 
damping ratios ("(tNi), are computed according to equations ( [3^ , with Nc = 5. 
Note that such damping ratios depend on the parameters m and s/m of the random 
field along with the hardening parameter H at meso-scale. For the sake of illus¬ 
tration, figure shows other results for another set of material parameters that is 
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Number of cycles Ni in free vibration 


Figure 14. Viscous-like damping ratio time history for Fg = 

5 kN (black) and Fq = 15 kN (grey). 

not optimal for representing the monotonic response in compression of the concrete 
used to build the tested column. The capability of the proposed material model 
for generating structural damping has been demonstrated and the development of 
an automatic procedure for identifying the full set of parameters targeting accurate 
representation of both cyclic concrete response and damping is left for future work. 




Figure 15. [top] Top-displacement time history in free vibration 
and [bottom] viscous-like damping ratio time history for a set of 
material parameters that is not optimal for the column considered 
above: m = 20 MPa, s/m = 10, C = 27.5 GPa and H = 10 GPa. 


6. Conclusions 

In this paper, a multi-scale stochastic uniaxial cyclic model suitable for rep¬ 
resenting most of the salient features of concrete nonlinear response observed in 
compressive experimental tests has been developed. It is based on the construction 
of a meso-scale where the response at each material point is elasto-plastic with 
kinematic hardening and heterogeneous yield stress. This implies that the transi¬ 
tion from elastic to plastic regime occurs at a loading level that is different in each 
material point. Heterogeneity is parameterized by a 2D homogeneous log-normal 
random field. As a first illustration of the capabilities of the model, some analyt¬ 
ical results are derived in the particular case of monotonic loading and vanishing 
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correlation length for the random field. Then, numerical simulations are performed 
and the effects of the parameters of the random field - that is the mean to, coef¬ 
ficient of variation s/to and correlation length - are investigated. It is shown 
that for small values of the correlation length, material response at macro-scale 
does not depend on the realization of the random field, showing that the devel¬ 
oped model is suitable for an objective representation of the material behavior. 
Besides, it is shown that the mean to and standard deviation s can be identified 
so that the monotonic compressive response of an actual concrete test specimen 
can be accurately represented by the developed model. The developed model how¬ 
ever lacks the ingredients for representing both strength and stiffness degradation 
mechanisms. Finally, the developed material model is implemented in a frame ele¬ 
ment in the purpose of representing the dynamic response of an actual reinforced 
concrete column. The numerical analysis of the column in free vibration shows the 
capability of the developed material model to create patterns classically associated 
to damping effects. In this simulation, damping does no come from some damping 
forces added in the dynamic balance equation (e.g. Rayleigh damping) but from 
the multi-scale stochastic nonlinear model. Although the underlying model is sto¬ 
chastic, the simulations and results shown are the same for any realization of the 
stochastic model. 

The main research prospects lie (i) in the enhancement of the model at meso- 
scale so that it can represent stiffness and strength degradation mechanisms at 
macro-scale; (ii) in the precise characterization of the stochastic model based on 
information from lower scales. This will consist in choosing, based on rational 
arguments, the type of first-order marginal law and correlation model, as well as 
the value of the corresponding parameters (mean, variance and correlation length). 
Although in another context, such an interaction between structural and material 
scientists has already been appealed for in [6] . Also, these issues could be considered 
in the context of stochastic micro-meso scale transition [431E1II2]. 
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